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ABSTRACT 

We report for the first time a 7 -ray and multiwavelength nearly-periodic oscillation in an active galactic 
nucleus. Using the Fermi Large Area Telescope (LAT) we have discovered an apparent quasi-periodicity in 
the 7 -ray fiux {E > 100 MeV) from the GeV/TeV BL Lac object PG 1553+113. The marginal significance of 
the 2.18 ± 0.08 year-period 7 -ray cycle is strengthened by correlated oscillations observed in radio and optical 
fiuxes, through data collected in the OVRO, Tuorla, KAIT, and CSS monitoring programs and Swift UVOT. The 
optical cycle appearing in ^ 10 years of data has a similar period, while the 15 GHz oscillation is less regular 
than seen in the other bands. Further long-term multi-wavelength monitoring of this blazar may discriminate 
among the possible explanations for this quasi-periodicity. 

Subject headings: gamma rays: galaxies — gamma rays: general — BL Lacertae objects: general — BL 
Lacertae objects: individual (PG 1553+113) — galaxies: jets — accretion, accretion disks 
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1. INTRODUCTION 

Among active galactic nuclei (AGN), blazars are distin¬ 
guished by erratic variability at all energies on a wide range 
of timescales. They are generally thought to be pow¬ 
ered by supermassive black holes (SMBHs, 10^-10^ Mq). 
P G 1553+113 GES 1553+113. z ^ 0.49 Danfo rth et 
al. I 201 QI : lAliu et al.l 120151 : 1 Abramowski et al.ll2015h is an 
optically/X-ray selected BL Lac object (iFalomo & TrevesI 
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1990 ) emitting variable GeV/T eV 7 radiation (lAleksic et al.l 
20151 : lAbramowski et al.ll2015h . As typical in very-high en¬ 
ergy (VHE) BL Lacs, the energetic non-thermal emission of 
PG 1553-fll3 originates in a relativistic jet and has a spec¬ 
tral energy distribution (SED) with two humps, overwhelm¬ 
ing any other component from either the nucleus or the host 
galaxy. 

The Large Area Telescope (LAT) on the Fermi Gamma-ray 
Space Telescope is providing continuous monitoring of the 
high-energy 7 -ray sky. The apparent modulation noted in the 
7 -ray flux of PG 1553-1-113 stimulated the multi-frequency 
and long-term variability study described in this paper. 

In 0 we describe the Fermi LAT data analysis and the 
sources of multiwavelength data; 0 details the multiple ap¬ 
proaches used for lightcurves and cross-correlation analysis; 
m outlines preliminary scenarios to interpret these results. 

2. FERMI LAT AND RADIO, OPTICAL, X-RAY DATA 

The LAT is a pair conversion detector with a 2.4 sr held of 
view, sens itive to y ravs from ^ 20 MeV to > 300 GeV l At- 
wood et a l.l2009l). The present work uses the new Pass 8 LAT 
database (I Atwood et al.ll2Ql^ . The LAT operating mode al¬ 
lows it to cover the entire sky every two ^ 1 . 6 -hour spacecraft 
orbits, providing a regular and uniform view of 7 -ray sources, 
sampling timescales from hours to years. This work uses ob¬ 
servations of PG 1553-i-113covering ^ 6.9 years (2008 Au¬ 
gust 4 to 2015 July 19, Modifled Julian Day, MJD, 54682.65- 
57222.65). The LAT data analysis employed the standard 
ScienceTools vlOrOpfQ package, selecting events from 
100 MeV-300 GeV with P8R2_SOURCE_V6 instrument re¬ 
sponse functions, in a circular Region of Interest of 10° ra¬ 
dius centered on the position of PG 1553-fll3. It used files 
gll_iem_v0 6 and iso_P8R2_SOURCE_V6_vO6 to model 
the Galactic and isotropic diffuse emission. Contamination 
due to the 7 -ray-bright Earth limb is avoided by excluding 
events with zenith angle > 90°. An unbinned maximum like¬ 
lihood model fit technique is applied to each time bin with a 
power-law spectral model and photon i ndex fixed to the 3F GL 
Catalog average value (1.604 ± 0.025, lAcero et ^12015h for 
PG 1553-1-113. The resulting lightcurves are shown in Fig. [T] 
Optical R-band data covering an interval of ^ 9.9 years 
(2005 April 19 to 2015 March 29, MJD 53479-57110) are 
reported in Fig. [21 Most unpublished observations were 
performed as part o f the Tuorla blazar monitoring program 
(iTakalo et al.l l2008 lf^. The data are reduced using a semi¬ 
automatic pipeline (Nilsson et al. in prep.). Public data from 
the Katzman Automatic Imaging Telescope (KAIT) and the 
Catalina Sky Survey (CSS) programs are also added. V-band 
magnitudes are scaled to the R-band values. 
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Fig. 1. — Fermi LAT 7 -ray lightcurves of PG 1553+113 over ~ 6.9 years, from 2008 August 4 to 2015 July 19. The lightcurve above 1 GeV is shown with a 
constant 45-day binning (top panel); two light curves above 100 MeV are shown, with 45- and 20-day binning (middle and bottom panels) 


As part of an ongoing blazar m onitoring program support¬ 
ing Fermi (iRichards et alJlMTTh . the Owens Valley Radio 
Observatory (OVRO) 40-m radio telescope has been observ¬ 
ing PG 1553+113 continually (about every 1 to 23 days) since 
2008 August. Figure [ 2 ] reports published 15 GHz lightcurves 
for the period from 2008 August 19 to 2014 May 18 (MJD 
54697-56795). OVRO in strumentation, data ca libration and 
reduction are described in IRichards et al.l (1201 Ih . 

Swift observed PG 1553+113 110 times between 2005 April 
20 and 2015 July 18 (unabsorbed 0.3-2 keV flux lightcurve in 
Figure [ 2 ]). X-Ray Telescope (XRT) data were first calibrated 
and cleaned (xrtpipeline, XRTDAS v.3.0.0) and energy 
spectra extracted from a region of 20 pixel (^47 ") radius, 
with a nearby 20 pixel radius region for background. Individ¬ 
ual XRT spectra are well fitted with a log-parabolic model, 
with c olumn density fixed to the Galactic value of 3.6 x 10^^ 
cm“^ (iKalberla et al.l201 Ih . Aperture photometry (5 " radius) 
for the UVOT V-band Alter was performed. 


3. TEMPORAL VARIABILITY ANALYSIS AND CROSS 
CORRELATION ANALYSIS 

We performed continuous wavelet transform (CWT) 
and Lomb-Scargle Periodogram (LSP) analyses on the 
lightcurves. Fig. [3] shows clear peaks at ^ 2 years for 7 - 
ray and optical power spectra. We also made an epoch folding 
(pulse shape) analysis used to extra ct the period, s hape, ampli¬ 
tude and phase, with uncertainties (iLarssonll 19961) . The for 
the folded pulse as a function of trial periods was fitted with 
a model containing 4 Fourier components, giving a period of 
798 ± 30 days (2.18 ± 0.08 years), consistent with the CWT 
and LSP findings (Fig. [3]). The value of the signal power peak 
does not change using r egular 20-day an d 45-day bins or an 
adaptive-bin technique (iLott et al.l 120121) for construction of 
the LAT lightcurve. 

A direct Power Density Spectrum (PDS) constructed from a 
LAT count-r ate lightcurve using exposure-w eighted aperture 
photometry (ICorbet et al.ll2007l : lKerill20lTI) above 100 MeV 
for a region with 3° radius with 600 second time bins (Fig. 
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Fig. 2.— Multifrequency lightcurves of PG 1553+113 at X-ray, optical and radio bands. Top panel: Swift-XRI integrated flux (0.3-2.0 keV). Central panel: 
optical flux density from Tuorla (R Alter, black fllled circle points), Catalina CSS (V filter rescaled, blue filled squared points), KAIT (V filter rescaled, red filled 
diamond points), Swift-XJWOT (V filter rescaled, green filled circle points). Dotted line: LAT flux (E > 100 MeV) with time bins of 20 days, scaled and y-shifted 
for comparison. Bottom panel: 15 GHz flux density from OVRO 40m (black filled circle points) and parsec-scale 15 GHz flux density by VLB A (MOJAVE 
program, yellow diamond filled points). 


©, confirms previous results with a peak at 2.16 ± 0.08 years, 
at 82 X the mean power level. The low-frequency modulation 
prevents an easy fit subtraction to the PDS continuum. The 
peak is ^ 5 times the mean level using a 4th order polynomial 
fit. 

The significance of any apparent periodic variation depends 
on what assumption is made about spurious stochastic vari¬ 
ability mimicking a periodic variation. The significance of 
the ^ 2 -year 7 -ray periodicity is difficult to assess given the 
limited length of the 7 -ray lightcurve. Red-noise, i.e. random 
and relatively enhanced low-frequency fluctuations over inter¬ 
vals comparable to the sample length, hinders the evaluation 
of per iodicity significance (e.g. iHsieh et al .120051 : ILaskv et all 
l2015h . We have approached the problem with two procedures: 

1) The red-noise is assumed to be produced by similar am¬ 
plitude flares (as seen in PG 1553-fll3 and some other LAT 
blazars), and the probability for these to line up in a regular 
pattern is estimated. The coherence of the periodic modula¬ 
tion was investigated by studying phase variations along the 
lightcurve. The local phase at each minimum and maximum 


was estimated by correlating a one-period long data segment 
with the Fourier template of the full lightcurve. The rms vari¬ 
ations relative to a perfectly coherent modulation was 27.4 
days. The chance probabilities for 3, 4 and 5 random events 
to be distributed with at least this coherence, as estimated by 
Monte Carlo simulations, are 0.0535, 0.0105 and 0.0027 re¬ 
spectively, implying a chance probability of a few percent for 
the 3.5-peak 7 -ray lightcurve of PG 1553-^113. 

2) We modeled the red-noise using Monte Carlo simu¬ 
lations with a first-order autoregressive process as the null 
hypothesis to assess whether the signal is consistent with a 
stochastic origin. Non-linear influence on the PDS is mini¬ 
mal thanks to the evenly spaced 7 -ray lightcurve. The power 
peak in Fig. [3]is above the 99% confidence contour level, i.e. 
has <1% chance of being a statistical fluctuation. The optical 
power peak has < 5% chance of being a statistical fluctuation. 

Although the 7 -ray periodicity signal alone is not com¬ 
pelling, the 9.9-years of optical data support the finding of 
a periodic oscillation in PG 1553-^113. The optical data, al¬ 
though affected by seasonal gaps, were analyzed using the 
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Fig. 3.— Left top panel: pulse shape (epoch-folded) 7 -ray (E > 100 MeV) flux lightcurve at the 2.18 year period (two cycles shown). Left bottom panels: 
2D plane contour plot of the CWT power spectrum (scalogram) of the 7 -ray lightcurve, using a Morlet mother function (fllled color contour). The side panel to 
this is the ID smoothed, all-epoch averaged, spectrum of the CWT scalogram showing a signal power peak in agreement with the 2.18-year value, also showing 
the LSP. Dashed lines depict increasing levels of confldence against red-noise calculated with Monte Carlo simulation. The 7 -ray signal peak is above the 99% 
confldence contour level (< 1% chance probability of being spurious). Right top panel: pulse shape from epoch folding of the optical flux lightcurve at the 2.18 
year period (two cycles shown). Right bottom panels: the same CWT and LSP diagrams for the optical lightcurve. The optical signal peak is above the 95% 
confldence contour level. 


same techniques as for the 7 -ray data. This analysis gives a 
period of 754 ± 20 days (2.06 ± 0.05 years), consistent within 
uncertainties with the 7 -ray results (Fig. El). 


f£ 


Period (days) 



Fig. 4.— Power Density Spectrum (PDS) of the LAT 0.1 — 300 GeV 
count rate lightcurve of PG 1553+113 from a 3° exposure-weighted aperture 
photometry technique with 600-second time bins. 


The less coherent 15 GHz lightcurve (5.7-years OVRO 
data) shows a signal power peak at 1.9 ± 0.1 year, with an 


additional power component at a 1.2-year timescale. Swift 
XRT data show a factor of 5 variation linearly correlated with 
the 7 -ray flux, while the synchrotron peak frequency shows a 
factor ^ 6 increase during high X-ray states, as suggested by 

iRdmeretalJd^ . 

The long-term X-ray count rate lightcurve from the Rossi- 
XTE ASM instrument (1996 February 20 to 2010 September 
11) and the Swift-BAI (from 2005 May 29) were also ana¬ 
lyzed but do not show any signal above the low-frequency 
noise, because of insufficient statistics. 

An important diagnostic for multi-frequency periodicity 
analysis is the discrete cross-correlation function (DCCF) 
used with two independent and complementary approaches. 

In the flrst procedure, flux variations are modeled assum¬ 
ing a simple power law oc 1//^ (with / = 1/t) in the PDS 
as measured directly from the lightcurve data, allowing us to 
estimate the cross-correlations signiflcance avoiding the as¬ 
sumption of equal variability in all sources a t the cost of a 
model assumption (iMax-Moerbeck et n]l2014l) . For the 7 -ray 
lightcurve with 20 -day binning we obtain a best lit o: = 0 . 8 , 
but the error is unconstrained, indicating that the length of 
the data set is too short (i.e. below flve cycles), relative to 
the suspected periodic modulation, to enable a reliable data 
characterization. The 45-day bin lightcurve yields a best lit 
(T = 0.1 with unconstrained error. The optical PSD is con¬ 
strained: the best fit value is a = 1.85, with la limits at 
[1.75, 2.00]. The 15 GHz flux light curve a slope of a = 1.4, 
with unconstrained limits on the a values as for the 7 -ray data. 
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Fig. 5.— Discrete cross-correlation plo ts from the approach with PDS 
model measured from the lightcurve data JMax-Moerbeck et "dll 20 1411 . In 
each plot the black dots are the DCCF estimates, and the red, orange and 
green lines are the la, 2a and 3a significance levels respectively. Top panel: 
DCCF between the radio 15GHz and 7 -ray (20-day time bins) lightcurve. 
Central panel: DCCF between the unbinned optical lightcurve and 7 -ray (20- 
day time bins) lightcurve. Bottom panel: DCCF between the 20-day rebinned 
optical lightcurve and 7 -ray (20-day time bins) lightcurve. The oscillating 
shape of the significance contours for this case is due to the number of sam¬ 
ples in each bin. 
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The DCCF between the unbinned radio lightcurve and the 20- 
day bin 7 -ray lightcurve results in a most probable time lag 
for radio-flux lagging the 7 -ray flux by 50 ± 20 days, with 
a 98.14% signiflcance for the best PSD flt with a range of 
[89.56%-99.99%] when flt erro rs are taken into account (Fig . 
[5]), using the fitting procedure of iMax-Moerbeck et al.l (120141) . 
The DCCF between the unbinned optical lightcurve and the 
20 -day bin 7 -ray lightcurve results in a most probable time 
lag for 7 -ray flux lagging the optical flux by 130 ± 14 days, 
with a 99.14% signiflcance for the best PSD flt and [96.09%- 
99.97%] when flt errors are taken into account (Fig. [5]). The 
DCCF peak is broad, however, and consistent with no lag. 
This is also seen when the optical data are rebinned into 20- 
day intervals, as shown in the bottom panel, where the most 
probable lag is 10 ± 51 days. 

In the second procedure, the signiflcance of the 7 -ray 
- radio correlation was estimated to be 95% by a mixed 


source correlation procedure (iFuhrmann et al.l [20141) . cross- 
correlating the PG 1553-fll3 lightcurve with those of 132 
comparison sources in that work, and evaluating the average 
DCCF level for time lags —100 to +100 days. The 7 -ray - 
optical correlation is significant at the 99% level, even though 
partly limited by the number of comparison sources and opti¬ 
cal lightcurve gaps. With only 132 comparison lightcurves we 
can measure a minimum probability-value of 0.0075, there¬ 
fore in principle a 99% level of signiflcance, but in this ap¬ 
proach the error in that estimate is hard to determine. With 
the mixed source methods there are two limitations: 1 ) the as¬ 
sumption that all the sources can be described with the same 
model for the variability, and 2 ) the sample variance due to the 
limited number of lightcurves must be assessed. The optical 
flux is found to lead the 7 -ray variations by 75 ± 27 days and 
the radio by 158 ± 10 days ( 7 -ray variations lead the 15 GHz 
flux variations by 83 ± 27 days). The possible reverse 7 -ray- 
optical time lag decreases to 28 ± 27 days when the optical 
lightcurve is binned. 

The possible opt ical- 7 -ray lag was already pointed out by 
ICohen et ^ (120141) . using KAIT unbinned optical lightcurves 
and LAT data. The high degree of 7 -ray-radio correlation 
in PG 1553-fll3 is not typically found in othe r individ¬ 
ual blazars/AGN (^see iMax-Moerbeck et akllMT^ . Signifi¬ 
cant cross-correlations are, neverthele ss, found when stackin g 
blazar samples (radio lagging 7 rays; IFuhrmann et al.|[ 20 T^ . 

4. DISCUSSION AND CONCLUSIONS 

Factors that led to the indication of a possible ~ 2-year 
periodic modulation in PG 1553-fll3 are: the continuous 
all-sky survey of Fermi’, the increased capability of the new 
Fermi LAT Pass 8 data; and the long-term radio/optical mon¬ 
itoring of 7 -ray blazars. Although the statistical signiflcance 
of periodicity is marginal in each band, the consistent positive 
cross-correlation between bands strengthens the case, mak¬ 
ing PG 1553-^113 the first possible quasi-periodic GeV 7 -ray 
blazar and a prime candidate for further studies. Hints of 
possible 7 -rfiy periodic ities are rare in literature (for example 
ISandrinelli et al.ll2014l) . The similarity of the low- and high- 
energy m odulation in PG 155 3-F113 is also a novel behavior 
for AGN (lRiegeill2()Q4l I2QQ7I) . Any periodic driving scenario 
should be related to the relativistic jet itself or to the process 
feeding the jet for this VHE BL Lac object. We outline, as 
examples, four possibilities: 

1. Pulsational accretion flow instabilities, approximating 
periodic behavior, are able to explain modulations in 
the energy outflow efficiency. Magnetically-arrested 
and magnetically-dominated accretion flows (MDAF) 
could be suitable regimes for rad iatively inefficient 
BL Lacs (iFragile & Meieil 120091) . characterized by 
advection-dominated accretion flows and subluminal, 
tu rbulent and peculiar radio kinematics (^ Karouzos et 
al. I2OI2I : iPiner & Edwards! 120141) . Su^ kinematics 
are sometimes explained a s a precessing or helical jet 
(IConwav & Murph^l 19931) . MDAF in a inner disc por¬ 
tion can be able to efficiently im part energy to parti¬ 
cles i n the jets of VHE BL Lacs (ITchekhovskov et ^ 
I2OIII) . Periodic instabilities are believe d to have short 
perio ds, ^ 10^ s • (M^mbi^/IO^ Mq) (iHonma et al.l 
1 19921) . but MHD simulations of magnetically choked 
accretion flows are se en to produce longer p eriods for 
slow-spinning SMBH (iMcKinnev et al.ll2QT^ . 

2. Jet precession (e.g., iRomero et al .1120001 : Stirling et 
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120041) . i.e. geometrical models 
presence of a jet wrapped by a sufficiently strong mag¬ 
netic field, could have a net apparent periodicity from 
the change of the viewing angle. Correspondingly the 
resulting Doppler magnification factor changes periodi¬ 
cally without the need for intrinsic variation in outfiows 
and efficiency. Non-ballistic hydrodynamical jet pre¬ 
cession m^explain variations with periods > 1 year 
(lRiegeill2004b . A differential Doppler factor AV{t) = 
r“^(l — P{t) cos6>(t))“^ < 40% variation (precession 
angle ^ 1 °) might be sufficient to support the ^ 2.8 
amplitude fiux modulation seen in 7 rays. A homoge¬ 
neous curve d helical jet scenari o for PG 1553+113 was 
proposed in iRaiteri et al.l (120151) . 


3. A mechanism analogous to low-frequency QPO from 
Galactic high-mass binaries/microquasars could pro¬ 
duce an accretion-outfl ow coupling mechanisrn as the 
basis of the periodicity (fpender & Bellonill2004l). King 
et al. (120131) ascribed the radio QPO in the FSRO 
CGRaBS J1359+4011 to this mechanism. However BL 
Lac objects like PG 1553+113 are thought to possess 
a lower accretion rate. The microq uasar QPO me cha- 
nism of Lense-Thirring precession (IWilkinslll972l) re¬ 
quires that the inner accretion flow forms a geometri¬ 
cally thick torus rather than a standard thin disc as the 
latter warps (^Bardeen-Petterson effect. Bardeen & Pet - 
terson ll975h rather than precesses (llngram et al.ll2009l) . 
A low mass accretion rate means that the accretion pro¬ 
cess probably forms an Advection-D ominated Accre¬ 
tion F low (ADAF), so it can process (iFragile & Meied 
l2009h . The X-ray emission in PG 1553+113 is prob¬ 
ably from the jet rather than from the flow, making it 
unlikely that the changing inclination of the hot flow 
causes the QPO. However, Lense-Thirring precession 
of the flow could affect the jet direction, giving the QPO 
as in ( 2 ) above. 


4. The pre sence of a gravitationally bound binary SMBH 
system (iBegelman et al.1 I198QI : iBarnes & Hernquls^ 
II 992 I) with a total mass ^10^ Mq, and a milli-pc sep¬ 
aration in the early inspiral gravitational-wave driven 
regime, might be another hypothesis. Keplerian binary 
orbital mo tion, would induce periodic accretion per¬ 
turbations (IValtonen et al.l 120081 : iPihajoki et al.l 120131 : 
iLiu et al.ll2015h or jet nutation expected from the mis¬ 
alignment of the rotating SMBH spins, or the gravi- 
tational torque on the disc exerted by the companion 
(lKatzlll997l:lRomero et all 120001 : ICaproni et all 120131 : 

iGraham et al.l2015(l . Significant acceleration of the disc 

evolution and accretio n onto a binary SMBH system is 
depic ted by modeling (iNixon et al.ll2013l : iDogan et^ 

Uml). 

Binary SMBH induced periodicities have timescales 


ranging from ^ 1 to ^ 25 years (iKomossal 120061 : 
lRiegeijl200^ . The SMBH total mass in PG 1553+113, 
estimated utilizing the putative link between in- 
fiow/accretion (d isc luminosity) and outfiow/jet (jet 
power) in blazars (iGhisellini et al.ll2014f) . is 1.6 x 10^ 
Mq, using a 0.1 MEdd rate and Doppler factor V = 30, 
in agreement with estimates for VHE BL Lacs (^ Woo et 
al. l2nn5h . 

The observed 2.18-year period is equivalent to an in¬ 
trinsic orbital time < Tohs/{^ z) 1.5 years, 
and the binary system size would be 0.005 pc (^ 100 
Schwarzschild radii). The probability to be observing 
such milli-pc system, estimated from the binary mass 
ratios 7 ^ 0.1 — 0.01 and the GW-driven regime lifetime 
(lPetersl[T964l) . tow — 10 ^ — 10 ^ years might be too 
small. 


Periodicities claimed for AGN are often controversial; 
however PG 1553+113 may potentially represent a key 7 - 
ray/multimessenger laboratory in the hypothesis of low- 
frequency gravitational wave e mission and may have asso - 
ciated PeV neutrino emission (iPadovani & RescomI 120141) . 
VLBI structure observations, radio/optical polarization data, 
and a prolonged multifrequency monitoring campaign will 
shed light on the situation. If the periodic modulation is real 
and coherent, as would be expected for a binary scenario, then 
subsequent maxima would be expected in 2017 and 2019, 
well within the possible lifetime of the Fermi mission. 
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